options(width = 255);
FILE_PATTERN = ".*cohort(20\\d{2})(current|lifespan)(\\d*)[^0-9].*"
FIELDS <- c("Article", "Article.Talk", "User", "User.Talk", "Other", "Maintenance");

# make sure to do this manually...
# install.packages(c("flexmix"), dependencies = TRUE)
library(flexmix);

main <- function() {
    lapply(c("dat/paper/brokendown", "dat/paper/brokendown-byuser"), function(dir) {
        lapply(c("hours", "edits", "bytes"), function(metric) {
            process_group(dir, metric);
        });
    });
};

process_group <- function(directory, metric) {
    paths <- Sys.glob(sprintf("%s/cohort*_prop_%s.stat", directory, metric));
    paths <- grep(FILE_PATTERN, paths, value = TRUE, perl = TRUE);
    write("=====================================================================", stdout());
    write(sprintf("processing %d files in %s with metric %s", length(paths), directory, metric), stdout());
    frames <- read_frames(paths);
    analyze_variation(frames);
    analyze_trends(frames);
};

read_frames <- function(paths) {
    return (lapply(paths, 
        function(path) {
            m <- regexec(FILE_PATTERN, path);
            parts <- regmatches(path, m);
            cohort <- parts[[1]][2];
            is_current <- parts[[1]][3];
            lifespan <- parts[[1]][4];
            #write(sprintf("processing file %s cohort %s, %s=%s", path, cohort, is_current, lifespan), stderr());
            data <- read.table(file=path, sep="\t", quote="", comment.char="", header=TRUE);
            data$cohort = as.integer(cohort);
            data$lifespan = as.integer(lifespan);
            data$is_current = (is_current == 'current');
            if (is_current == 'current') {
                data$lifespan = "c";
            }
            data <- data[data$Date != "2012",];
            return (data);
        }
    ));
};


analyze_trends <- function(frames) {
    all_data <- Reduce(function(x, y) merge(x, y, all=T), frames, accumulate=F)

    z <- lapply(FIELDS, function(f) {
        m <- aov(all_data[,f] ~ factor(all_data$cohort) + factor(all_data$lifespan) + factor(all_data$Date));
        print(f);
        print(summary(m));
        print(drop1(m));
        print("");
        print("");
        print("");
    });

    #Y = cbind(all_data$Article, all_data$Article.Talk, all_data$User, all_data$User.Talk, all_data$Other);
    #fit <- manova(Y ~ factor(all_data$cohort) + factor(all_data$lifespan) + factor(all_data$Date));
    #print(fit);
    #print(summary(fit, test="Pillai"));
    #print(summary.aov(fit));
}

analyze_variation <- function(frames) {
    frames <- Filter(function(f) nrow(f) > 1, frames);
    kls <- sapply(frames, function(f) {
        f <- f[,FIELDS];
        kl <- KLdiv(t(f));
        n <- nrow(f);
        m <- sum(kl) / (n*n - n);   # mean except fo diagonals (which are 0)
        return (m);
    });
    
    write(sprintf("kl_div n=%d mean=%f dev=%f", length(kls), mean(kls), sd(kls)), stdout());

    lapply(FIELDS, function(field) {
        means <- sapply(frames, function(frame) { return (mean(frame[field])); });
        sds <- sapply(frames, function(frame) { return (sd(frame[field])); });
        write(sprintf("ns %s n=%d mean(mean)=%f mean(sd)=%f",
                field, length(sds), mean(means), mean(sds)), stdout());
    });
}

foo <- main();
